Naval  Research  Laboratory 


ft  M 


Washington,  DC  20375-5320 


NRL/MR/7643--95-7772 


Calculation  of  the  Ionospheric  0+ 
Concentration  From  O II 834  A  Airglow 
Using  Discrete  Inverse  Theory 


J. M.  PlCONE 
R.R.  Meier 

K. F.  Dymond 
R.P.  McCoy 

E.O.  Hulburt  Center  for  Space  Research 
O.  Kelley 

George  Mason  University 
Faiifax,  VA  22030 

R.J.  Thomas 

New  Mexico  Institute  of  Mining  and  Technology 
Socorro,  NM  87801 


September  1,  1995 


OTIC 

ELECTEE 

SEP  2  6  1995 1 


19950922  124 


Approved  for  public  release;  distribution  unlimited. 


Form  Approved 
OMB  No.  0704-0188 


REPORT  DOCUMENTATION  PAGE 


Public  reporting  burden  for  this  collection  of  informetion  is  estimeted  to  everege  1  hour  per  response,  including  the  time  for  reviewing  inetrui 
gathering  end  maintaining  the  date  needed,  end  completing  end  reviewing  the  collection  of  informetion.  Send  comments  regarding  this  burde 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  C 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduction  Project  (0704 

stiorw,  Marching  existing  data  sources, 
n  estimate  or  any  other  aspect  of  this 
)perations  and  Reports,  1216  Jefferson 
-0188),  Washington,  DC  20603. 

1.  AGENCY  USE  ONLY  {Leave  Blank) 

2.  REPORT  DATE 

September  1,  1995 

3.  REPORT  TYPE  AND  DATES  COVERED 

Memo  Report  1  Oct  1991-1  July  1995 

4.  TITLE  AND  SUBTITLE 

Calculation  of  the  Ionospheric  0+  Concentration  From  0  II  834  A  Airglow  Using  Discrete 

Inverse  Theory 

5.  FUNDING  NUMBERS 

PE-61 153N 

PE-63716D 

PE-35 160F 

6.  AUTHOR  (SI 

J.M.  Picone,  R.R.  Meier,  O.  Kelley*  K.F.  Dymond,  R.J.  Thomas,**  and  R.P.  McCoy 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Research  Laboratory 

Washington,  DC  20375-5320 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

NRL/MR/7643 -95-7772 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Office  of  Naval  Research 

800  North  Quincy  Street 

Arlington,  VA  22217-5660 

10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES 

♦Institute  for  Computational  Sciences  and  Informatics,  George  Mason  Univ.,  Fairfax,  VA  22030 
♦♦Department  of  Electrical  Engineering,  New  Mexico  Institute  of  Mining  and  Tech.,  Socorro,  NM  87801 

12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited. 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  {Maximum  200  words ) 

Discrete  inverse  theory  (DIT)  forms  the  basis  of  new  techniques  for  extracting  dayside  0+  number  density  profiles  from  the  834  A 
airglow,  as  measured  by  a  limb-scanning  system  on  an  orbiting  satellite.  Our  tests  of  this  method  assume  observations  from  an 

altitude  of  850  km  with  scans  from  10°  to  26.5°  below  horizontal,  consistent  with  future  multiyear  missions.  The  retrieval  code 
computes  an  iterative,  maximum  likelihood  solution  by  comparing  observations  to  estimates  calculated  with  a  new  forward  model. 
The  model  includes  multiple  resonant  scattering  and  pure  absorption.  To  generate  synthetic  data  for  tests,  we  represent  the  "true"  0+ 
distribution  as  a  Chapman  layer  and  compute  an  intensity  profile  wdth  the  forward  model,  adding  simulated  noise.  We  present 
detailed  studies  of  convergence  properties  of  the  retrieval  techniques  and  of  uncertainties  in  the  retrieved  parameters.  For  this 
baseline  (Chapman  layer)  case,  the  method  is  robust,  converging  to  an  accurate  solution  for  a  wide  variation  in  synthetic  data.  We 
include  a  brief  preview  of  recent  studies  showing  the  following:  (1)  for  future  missions,  the  DIT  method  can  correctly  distinguish 
between  distinctly  different  Chapman  layers  that  produce  nearly  identical  intensity  profiles  and  (2)  the  retrieval  of  an  additional 
parameter  that  scales  the  model  intensity  profile  can  compensate  for  inaccuracies  in  the  instrument  sensitivity  or  in  the  magnitude  of 
the  initial  volume  excitation  rate. 

14.  SUBJECT  TERMS 

Remote  sensing  Discrete  inverse  theory 

Computer  algorithms  Upper  atmosphere 

Ionosphere 

15.  NUMBER  OF  PAGES 

31 

16.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION 

OF  REPORT 

UNCLASSIFIED 

18.  SECURITY  CLASSIFICATION 

OF  THIS  PAGE 

UNCLASSIFIED 

19.  SECURITY  CLASSIFICATION 

OF  ABSTRACT 

UNCLASSIFIED 

20.  LIMITATION  OF  ABSTRACT 

UL 

NSN  7540-01-280-5500  Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std  239-18 


298-102 


CONTENTS 


1.  Introduction . 

2.  Inversion  Technique . 

2.1.  Discrete  Inverse  Theory  (DIT) . 

2.2.  Estimating  Uncertainties  in  Retrieved  Model  Parameters 

2.3.  Estimation  of  O II  834  A  Limb  Intensities . 

3.  Model  System  and  Simulated  Observations  for  Tests.... 


4.  Retrieval  of  Chapman  Parameters . 

4. 1 .  Typical  Results  for  the  Three-Parameter  Model . 

4.2.  Convergence  Properties . 

4.3.  Uncertainties  in  Retrieved  Model  Parameters . 

4.3.1.  Uncertainty  in  Retrieved  Nmax . 

4.3.2.  Uncertainty  in  Retrieved  zmax . 

4.3.3.  Uncertainty  in  Retrieved  Scale  Height . 

4.4.  Uncertainties  in  Other  Components  of  the  Forward  Model 


12 

12 

13 

.14 

15 

16 
.16 
.17 


5.  Concentration  Scalar  Model 


18 


6.  Summary  and  Discussion:  Practical  Applications . 

6.1.  Assumptions  Regarding  the  Observations  and  the  Ionosphere . 

6.2.  Numerical  Procedures  and  Solution  Criteria . 

6.3.  Advanced  Topics . 

6.3.1.  Different  Chapman  Layers  With  "Nearly  Identical"  Intensity  Profiles 

6.3.2.  Retrieval  of  a  Data  Scaling  Factor . 

7.  Acknowledgments . 

Appendix:  Numerical  Details  of  the  Forward  Model . 

References . . . 


.20 

.21 

.21 

.22 

.23 

.24 

.25 

.26 

.28 


Figures 


,31 


Accesion  For 

/ 

NTIS  CRA&I 

DTIC  TAB 

□ 

Unannounced 

Justification 

□ 

By . . . 

Distribution/ 

Availability 

Codes 

Aval!  and  /  or 
Special 


CALCULATION  OF  THE  IONOSPHERIC  O+  CONCENTRATION  FROM 
O  n  834  A  AIRGLOW  USING  DISCRETE  INVERSE  THEORY 


1.  Introduction 

The  0+  (2s^  2p^  ^S°  <—  2s  2p^  ^P)  transition  produces  a  triplet  emission  at  line  center 
wavelengths  of  832.76,  833.33,  and  834.47  A  [Eriksson,  1987].  This  is  the  brightest  O  H  feature 
observed  in  the  ultraviolet  (UV)  dayglow  [Meier,  1991].  The  834  A  triplet  transition  is  an 
allowed  multiplet,  for  which  the  photons  undergo  multiple  resonant  scattering  by  0+  ions  in  the  F 
layer  of  the  ionosphere.  This  process  significantly  enhances  the  834  A  volume  emission  rate  in 
the  vicinity  of  the  F-layer,  so  that  different  0+  distributions  can  produce  qualitatively  different  F- 
layer  limb  intensity  profiles  [McCoy  et  al.,  1985].  The  purpose  of  the  present  work  is  to 
demonstrate  new  techniques  for  the  retrieval  of  altitude  profiles  of  the  0+  number  density  from 
observations  of  the  O  H  834  A  airglow  by  a  limb-scanning  extreme  ultraviolet  (EUV) 
spectrograph.  Because  0+  is  the  dominant  ion  in  the  F  region  of  the  ionosphere,  these  methods 
can  transform  834  A  limb  intensity  data  into  electron  density  maps  of  the  F  region. 

This  work  supports  a  number  of  multiyear  satellite-borne  systems  for  remote  sensing  of 
the  upper  atmosphere;  the  observations  will  take  place  over  the  next  fifteen  years  and  should 
continue  well  beyond  that  period.  These  systems  include  the  Remote  Atmospheric  and 
Ionospheric  Detection  System  (RAIDS),  for  which  a  new  host  platform  has  not  yet  been 
identified  [McCoy  et  al.,  1994a;  Christensen  et  al.,  1993];  the  Advanced  Research  and  Global 
Observation  Satellite  (ARGOS;  McCoy  [1995]),  due  for  launch  in  1996,  and  the  Special  Sensor 
Ultraviolet  Limb  Imager  (SSULI:  McCoy  et  al.  [1994b]),  which  will  be  flown  by  the  Defense 
Meteorological  Satellite  Program  (DMSP)  on  a  series  of  five  missions  possibly  beginning  as 
early  as  1997.  All  of  these  systems  will  occupy  polar,  Sun-synchronous,  circular  orbits  at 
altitudes  between  800  and  900  km  and  will  scan  the  limb  of  the  Earth  within  an  angular  range  of 
approximately  10°  to  26.5°  below  the  horizontal.  All  missions  carry  instruments  that  provide 
spectral  resolution  of  12  to  15  A  in  the  EUV  band  and  7  to  1 1  A  in  the  FUV  band.  Given  the 
similarity  of  the  respective  measurements,  these  missions  will  provide  a  long-term,  self- 
consistent  data  set  for  detailed  global  studies  of  the  dynamical  behavior  of  the  upper  atmosphere 
as  a  function  of  various  geophysical  parameters  and  forcing  agents,  e.g.,  solar  and  geomagnetic 
activity.  In  addition,  the  database  should  be  useful  in  developing  more  accurate  empirical  and  ab 
initio  models  of  the  upper  atmosphere. 

McCoy  et  al.  [1985]  laid  the  foundation  for  these  developments  by  performing  detailed 
parametric  studies  of  the  relationships  between  various  geophysical  and  ionospheric  parameters 
and  O II  834  A  airglow  measurements  from  a  spinning  satellite  (STP78-1).  The  work  relied 
primarily  on  a  radiative  transport  model  that  included  the  effects  of  multiple  resonant  scattering 
by  o+  and  pure  absorption  (extinction)  of  photons  by  O,  N2,  <md  O2  [Anderson  and  Meier, 
1985].  To  determine  the  location  and  rate  of  photon  production,  including  both  photoionization 
and  photoelectron  impact  ionization  of  atomic  oxygen,  McCoy  et  al.  used  the  PEGFAC  code  of 
Strickland  and  Meier  [1982],  A  standard  Chapman  layer  represented  the  0+  distribution  in  the 
altitude  variable  z : 
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where  NQ+  is  the  0+  number  density  at  a  given  altitude,  zmax  is  the  altitude  at  which  the  0+ 

number  density  peaks,  Nmax  is  the  peak  0+  number  density,  and  H  is  the  (neutral)  atomic  oxygen 
scale  height.  The  study  indicated  that  the  solar  flux,  the  photoelectron  flux,  and  the  neutral 
number  density  distributions  mainly  affected  the  magnitude  of  the  intensity  profiles  while  the 
shape  changed  little.  On  the  other  hand,  changes  in  the  Chapman  parameters  [Nmax  ,  zmax  ,  H] 
caused  distinct  changes  in  the  shape  of  the  intensity  profile,  suggesting  that  these  parameters 
could  be  retrieved  from  the  data. 


The  present  paper  presents  a  thorough  analysis  of  efficient  methods  for  retrieving  an 
optimal  0+  number  density  distribution  from  834  A  limb  intensity  measurements.  The  basis  for 
the  methods  is  discrete  inverse  theory  (DIT:  Rodgers  [1976],  Tarantola  [1987],  Tarantola  and 
Valette  [1982],  Menke  [1989]),  which  has  been  introduced  by  Meier  and  Picone  [1994]  for  the 
retrieval  of  thermospheric  composition  and  temperature  from  far  ultraviolet  (FUV)  limb  intensity 
profiles.  To  establish  a  direct  connection  to  the  work  of  McCoy  et  al.  [1985]  and  to  provide  an 
assessment  of  the  potential  of  DIT  for  moderately  optically  thick  systems,  we  have  used  the 
Chapman  layer  representation  of  the  0+  distribution  in  the  F  region,  as  given  by  Equation  (1). 
This  representation  has  been  useful  for  both  qualitative  and  quantitative  analysis  of  the  F  region 
[Kelley  and  Heelis,  1989;  Chamberlain  and  Hunten,  1987],  On  the  other  hand,  Anderson  et  al. 
[1987]  have  pointed  out  that  Equation  (1)  is  not  general  enough  to  provide  close  fits  to  profiles 
produced  by  sophisticated  simulation  models.  D.  N.  Anderson  and  M.  Fox  [private 
communication,  1 993]  have  suggested  generalizations  of  the  Chapman  layer  with  additional 
parameters.  The  present  paper  will  provide  a  useful  baseline  for  assessing  more  robust 
representations  of  the  0+  distribution. 

The  specific  task  covered  by  this  paper,  therefore,  is  to  retrieve  the  three  Chapman 
parameters  underlying  synthetic  limb  intensity  data  that  are  consistent  with  the  sensitivity  and 
observational  geometry  of  the  systems  mentioned  above.  We  follow  the  presentation  of  Meier 
and  Picone  [1994]  on  the  retrieval  of  thermospheric  composition  and  temperature  from  the  FUV 
dayglow.  The  key  difference  is  that  Meier  and  Picone  neglected  resonant  scattering  effects, 
while  such  effects  are  significant  in  the  case  of  the  O II  834  A  field.  Multiple  scattering  can 
enhance  the  volume  emission  rate  in  the  F-layer;  on  the  other  hand,  a  sizable  loss  of  information 
can  occur  along  the  line  of  sight.  These  effects  increase  the  complexity  and  nonlinearity  of  the 
forward  model. 


Section  2  presents  a  summary  of  discrete  inverse  theory  and  discusses  the  forward  model 
for  calculating  the  limb  intensity  during  the  retrieval.  We  use  the  same  model  to  generate 
synthetic  data  as  well.  The  Appendix  outlines  the  numerical  algorithm  used  for  computing 
volume  excitation  rates  according  to  the  prescription  in  Section  2.  Section  3  describes  the  model 
system  and  the  simulated  observations  used  for  testing  the  DIT  code.  The  tests  employ  the 
representation  of  the  0+  number  density  given  by  Equation  (1)  both  to  generate  synthetic  data 
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and,  as  part  of  the  DIT  forward  model,  to  estimate  intensity  values  during  the  retrieval  process. 
As  a  result,  the  tests  assume  that  the  DIT  forward  model  constitutes  an  accurate  representation  of 
the  true  physical  system.  Consequently,  the  results  provide  an  upper  bound  on  the  performance 
of  DIT  for  this  application. 


Section  4  discusses  the  convergence  properties  of  the  algorithms  and  the  sensitivity  to  the 
initial  conditions  for  the  iterative  process.  We  then  provide  a  picture  of  the  uncertainties  in  the 
retrieved  Chapman  parameters  as  functions  of  the  "true"  parameter  values  used  to  generate  the 
synthetic  data.  Section  5  uses  the  "concentration  scalar  model"  (CSM)  of  Meier  and  Picone 
[  1 994]  and  R.  Thomas  [private  communication,  1 993]  to  assess  the  sensitivity  of  limb-scanning 
experiments  to  the  0+  density  at  different  altitudes  within  the  F-layer.  Section  6  summarizes  the 
main  results  of  the  study  and  discusses  the  details  involved  in  producing  comprehensive, 
practical  O II  834  A  inversion  algorithms  which  are  customized  to  various  scientific  and 
engineering  applications.  Because  of  the  interest  generated  by  recent  conference  presentations 
[Picone  et  al.,  1994;  Kelley  et  al.,  1994]  and  publications  [Link  et  al.,  1994],  Section  6  includes 
a  brief  preview  of  ongoing  studies  [Picone  et  al.,  1995a,b]  addressing  the  following  issues. 

(1)  the  retrieval  of  an  additional  parameter  to  scale  the  intensity  data  and 

(2)  the  problem  of  faithfully  retrieving  distinctly  different  0+  number  density  profiles  which 
correspond  to  "nearly  identical"  intensity  profiles. 
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2.  Inversion  Technique 


2.1.  Discrete  Inverse  Theory  (DIT) 

Meier  and  Picone  [1994]  have  reviewed  the  primary  features  of  Menke’s  [1989] 
presentation  of  DIT  with  nonlinear  forward  models.  This  formulation,  which  is  based  on  a 
maximum  likelihood  criterion,  applies  specifically  to  cases  in  which  the  distributions  of  data 
points  and  model  parameters  about  their  respective  mean  values  are  approximately  Gaussian. 
Here  we  briefly  describe  and  extend  the  discussion  of  Meier  and  Picone  to  include  the 
relationship  of  Menke’s  approach  to  others. 

Consider  a  set  of  discrete  observational  data  which  may  be  expressed  as  components  of  a 
vector  d°  of  length  N.  In  the  present  context,  this  vector  contains  limb  intensity  measurements  as 
a  function  of  observation  zenith  angle  (OZA).  In  order  to  extract  information  regarding  the 
concentration  of  the  emitting  species,  we  require  a  forward  model, 

G(m)  =  d,  (2) 

in  which  m  is  a  vector  of  parameters  defining  the  spatial  distribution  of  the  emitting  species  and 
the  vector  d  is  the  "model"  limb  intensity,  computed  as  a  function  of  OZA  for  a  given  m.  The 
vector  function  G  thus  maps  a  species  density  distribution  onto  a  limb  intensity  profile 
corresponding  to  a  particular  scan.  We  denote  the  number  of  parameters,  and  hence  the  length  of 
m,  by  the  variable  M.  For  the  present  tests  of  our  algorithms,  the  number  of  parameters,  M,  is 
three  and  m  =  [Nmax  ,  zmax  ,  H],  the  vector  containing  the  Chapman  parameters  in  Equation  (1). 

As  we  will  show  in  the  next  subsection,  the  model  limb  intensity  vector  G(m)  is  a  highly 
nonlinear  function  of  the  0+  number  density  and  the  major  neutral  constituents.  To  compute 
parameter  values,  we,  therefore,  use  an  iterative  technique  equivalent  to  that  of  Meier  and  Picone 
[1994].  We  seek  the  maximum  likelihood  solution  for  cases  in  which  both  the  data  and  model 
parameters  are  distributed  normally,  i.  e., 

P(d,d°)  -  exp  {-  }  (d-d°)T  [cov  d]  -1  (d-d0)}  (3) 

and 

P(m,(m))  °c  exp  ^  (m-(m})T  [covm]-1  (m-<m>)}  .  (4) 

In  Equation  (3),  [cov  d]  is  the  covariance  matrix  for  the  data.  In  Equation  (4),  (m)  is  the  vector 
of  a  priori  model  parameter  values,  which  provide  estimates  of  the  parameter  values  to  be 
retrieved,  and  [cov  m]  is  the  covariance  matrix  of  the  a  priori  model  parameter  values.  The  a 
priori  model  parameter  vector  permits  the  calculation  of  the  optimum  estimate  of  model 
parameters  when  the  problem  is  underdetermined  or  mixed-determined.  Possible  sources  for 
values  of  (m)  include  results  from  the  analysis  of  prior  data  sets  (e.g.,  preceding  limb  scans), 
semi-empirical  models,  or  theoretical  constraints  or  estimates. 
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To  obtain  the  maximum  likelihood  solution  for  a  given  set  of  observations  and  for  the 
corresponding  a  priori  model  parameter  values,  we  compute  the  best  estimate  me  of  the  vector  m 
which  maximizes  the  joint  probability 

P(m,  d)  oc  P(m,  <m»  P(d,  d°)  ,  (5) 

subject  to  the  constraint  of  Equation  (2).  Hence  we  must  minimize  the  generalized  %2  >  given  by 

X2  =  (d-d°)T  [cov  d] -1  (d  -  d° )  +  (m  -  (m))T  [cov  m] _1  (m  -  (m))  .  (6) 

One  well-known  method  of  accomplishing  this  is  the  Levenberg-Marquardt  method  (LM:  e.g., 
see  Press  et  al.  [1989]).  Meier  and  Picone  [1994]  presented  an  alternative  method  based  on  the 
formulation  of  Menke  [1989]  (also  Tarantola  [1987]).  The  two  methods  provide  a  solution  when 
the  forward  model,  Equation  (2),  is  nonlinear.  Both  start  with  an  initial  estimate,  me0,  and 
iteratively  solve  an  approximate  equation  in  which  the  solution  at  step  (n  +  1)  depends  both  on 
men  ,  the  estimate  of  the  model  parameters  after  the  nth  iteration,  and  on  the  gradient,  VGn  ,  of 
the  forward  model  with  respect  to  m,  evaluated  at  m  =  men.  Thus,  no  higher  derivatives  of  the 
forward  model  occur  in  the  procedures.  Unlike  the  calculations  of  McCoy  et  al.  [1985],  the 
methods  allow  all  parameters  to  change  at  each  computational  step.  In  our  implementation,  the 
only  numerical  approximation  of  note  is  the  use  of  one-sided  derivatives  to  estimate  the  gradient 
matrix  VGn  .  Higher-order  derivatives  require  that  the  forward  model  be  evaluated  at  more 
points  in  parameter  space.  The  concomitant  increased  expense  has  proven  to  be  more  significant 
than  any  loss  of  accuracy  in  the  search  for  an  optimum  solution.  For  the  tests  in  this  paper,  the 
procedure  has  "converged"  when  the  change  in  successive  values  of  x2  is  10'4  or  less. 

Comparison  with  results  based  on  more  stringent  convergence  criteria  have  shown  that  the  use  of 
10‘4  translates  into  an  uncertainty  in  the  range  of  ±  2%  for  me . 

From  the  presentation  of  LM  by  Press  et  al.  and  the  discussion  of  the  alternative  method 
by  Tarantola  [1987],  we  have  found,  in  fact,  that  the  two  formulations  (i.e.,  Menke  and  LM)  are 
equivalent  when  applied  to  the  generalized  x2  of  Equation  (6).  (Note  that  Press  et  al.  deal 
specifically  with  the  overdetermined  case  and  do  not  allow  for  a  priori  data  with  a  finite 
covariance.)  Tarantola  calls  this  a  "quasi-Newton  method,"  which  has  superior  convergence 
properties  in  comparison  to  the  method  of  steepest  descent.  The  only  formal  difference  between 
the  two  approaches  lies  in  selecting  the  size  of  the  change  in  the  model  parameter  vector  from 
iteration  to  iteration.  As  the  optimum  value  of  me  is  approached,  one  would  prefer  to  take 
smaller  steps,  because  of  the  numerical  and  analytic  approximations  that  are  involved.  If  the  step 
size  is  too  large  on  the  nth  iteration,  for  example,  the  nth  estimate,  men  ,  might  actually  be  farther 
from  the  region  of  convergence  than  is  men.i-  Thus  satisfactory  convergence  might  never  be 
achieved.  LM  checks  for  this  specific  situation,  adjusting  the  step  size  iteratively  to  ensure  better 
convergence  at  each  step.  Menke  includes  no  corresponding  calculation  of  optimal  step  size. 
Fortunately  the  user  can  incorporate  the  LM  step  size  selection  procedure  into  Menke’s  equations. 
Tarantola  also  provides  analytic  methods  of  choosing  an  appropriate  step  size;  we  have  not 
investigated  his  analytic  estimates. 
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2.2.  Estimating  Uncertainties  in  Retrieved  Model  Parameters 


Of  nearly  equal  importance  with  the  retrieval  of  model  parameters  is  the  ability  to 
estimate  the  uncertainties  in  both  the  retrieved  parameter  values  and  the  values  of  the  physical 
quantities  which  have  been  parameterized.  In  our  case,  the  parameterized  physical  quantity  is  the 
0+  number  density  as  a  function  of  altitude.  Here  Menke  [1989]  and  Tarantola  and  Valette 
[1982]  transcend  the  LM  approach  by  providing  an  estimate  of  the  covariance  matrix  for  me, 
based  on  the  results  for  the  general,  linear  Gaussian  case: 

[coving]  =  VGjfv[cov  d]  (VGjfv  )T  +  [l-Rn]  [cov  m]  [i  —  R n ] T 

,  x  ,  A  ,  •  (?) 

=  (VGnT[covd]  1  VGn  +  [covm]  ) 

In  Equation  (7),  men  is  the  estimate  for  the  optimum  model  parameter  vector  after  the  n^1 
iteration,  VGinvn  is  the  "generalized  inverse  matrix"  for  the  gradient  VGn  of  the  forward  model, 
I  is  the  identity  matrix  for  model  parameter  space,  and  the  "model  resolution  matrix"  is 
[Tarantola,  1987;  Menke,  1989]: 


Rn  -  VGf  VGn  .  (8) 

For  an  overdetermined  case  (e.g.,  Chapman  layer  representation  with  N  >3),  Rn  =  I .  On  the 
other  hand,  the  retrieval  problem  is  often  mixed-determined  [Menke,  1989]  when  the 
concentration  scalar  model  is  used  (Section  5),  and  Rn  *  I .  The  standard  deviation  Oj  = 

([cov  me]n)'/2  gives  a  measure  of  the  uncertainty  in  the  retrieved  model  parameter  mej . 

If  the  forward  model  equations  or  their  numerical  solutions  are  approximate,  or  if  the 
forward  model  depends  on  quantities  that  are  not  known  precisely  (e.g.,  the  values  of  cross 
sections),  Menke  [1989]  and  Tarantola  [1987]  represent  the  model  estimates  of  d  by  a  Gaussian 
distribution  about  G(m)  with  a  covariance  matrix,  [cov  G],  that  corresponds  to  the  uncertainty  in 
the  estimates.  Equation  (5)  then  must  include  this  distribution  as  a  third  multiplicative  factor, 
rather  than  assuming  the  exact  constraint  of  Equation  (2).  Equation  (6),  Equation  (7),  and  the 
definition  of  VGmvn  then  hold  with  [cov  d]  replaced  by  [cov  d]  +  [cov  G]  [Tarantola  and 
Valette,  1982] .  These  modifications  permit  an  assessment  of  the  effects  of  approximations  in 
the  forward  model  and  of  uncertainties  in  inputs  to  the  forward  model  on  the  overall  uncertainty 
in  the  retrieved  model  parameters.  The  intent  of  the  present  paper  is  to  provide  a  baseline  for 
such  assessments;  thus  we  assume  that  the  constraint  of  Equation  (2)  is  exact. 

To  convert  the  uncertainties  in  the  model  parameters,  as  represented  by  [cov  me],  to  the 
uncertainty  in  the  0+  number  density  at  a  given  altitude,  we  use  the  standard  formula  for 
propagation  of  errors  [Bevington  and  Robinson,  1992]: 
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var[No+(z) 


(9) 


"3No+(z)^ 

f<aNo+(z)>| 

3mi  J 

l  amJ  J 

Here  the  term  "var"  represents  the  variance,  i.e.,  the  square  of  the  standard  deviation.  We 
evaluate  the  derivatives  using  Equation  (1)  with  the  retrieved  model  parameter  values. 

Since  the  covariance  matrix  is  symmetric,  [cov  me]  has  ViiM1  +  M)  independent 
elements.  Equation  (9),  therefore,  can  provide  realistic  estimates  of  the  uncertainties  at  each  of  a 
large  number  of  grid  points  only  if  Equation  (1)  accurately  represents  the  true  0+  number  density 
distribution.  In  that  case,  Equation  (1)  can  provide  a  faithful  extrapolation  of  the  0+  number 
density  to  altitudes  which  are  not  well-sampled  by  the  observing  instrument. 

To  determine  the  amount  of  information  that  the  instrument  is  receiving  from  each  region 
of  the  upper  atmosphere,  we  may  generate  synthetic  data  using  a  Chapman  layer  or  another 
model  and  then  use  the  concentration  scalar  model  (CSM)  to  represent  the  0+  number  density 
profile  in  the  retrieval  code  [Meier  and  Picone,  1994].  The  CSM  specifies  a  set  of  M(<N) 
altitude  grid  points  at  each  of  which  the  0+  number  density  is  itself  treated  as  a  model  parameter 
independent  of  the  number  density  at  any  other  altitude  point.  In  contrast  to  Equation  (1)  or 
semi-empirical  models,  the  densities  at  different  grid  points  are  not  coupled.  The  CSM  number 
density  profile,  therefore,  cannot  extrapolate  values  from  the  region  actually  providing 
observational  information  to  other  regions  which  are  not.  This  modeling  approach,  described  in 
Section  5,  will  reveal  the  degree  of  uncertainty  in  the  0+  number  density  at  a  given  altitude, 
arising  from  the  loss  of  information  that  is  being  transferred  along  the  line  of  sight. 

2.3.  Estimation  of  O II  834  A  Limb  Intensities 


The  0+  ion  is  the  dominant  ion  in  the  F  layer  of  the  ionosphere,  and  assuming  charge 
neutrality,  knowledge  of  the  0+  distribution  there  constitutes  knowledge  of  the  electron  number 
density  there  as  well.  Meier  [1991]  has  discussed  the  contribution  of  O II 834  A  emissions  to  the 
airglow  in  the  upper  atmosphere  and  the  relationship  of  0+  to  the  ionosphere  and  magnetosphere, 
including  experimental  measurements.  According  to  that  assessment,  834  A  emissions  can  be 
moderately  optically  thick,  primarily  due  to  resonant  scattering  by  0+  of  photons  originally 
produced  below  200  km  by  photoionization  and  photoelectron  impact  ionization  of  atomic 
oxygen.  Link  et  al.  [1994]  have  discussed  laboratory  measurements  of  cross  sections  for  O II  834 
A  production,  resonant  scattering,  and  absorption  and  have  related  the  results  to  the  airglow.  We 
have  incorporated  those  results  into  our  forward  model  for  predicting  limb  intensities.  Under  the 
conditions  used  in  our  tests  (Section  3)  and  for  a  Chapman  layer  with  [Nmax  ,  zmax  ,  H]  =  [1.0  x 
106  cm-3,  350  km,  60  km],  the  maximum  vertical  optical  depth  for  resonant  scattering  of  the 
834.47  A  line  is  approximately  4.2  (z  *  216  km).  The  corresponding  vertical  optical  depth  for 
pure  absorption  (photon  loss  through  absorption  by  N2, 02,  and  O)  is  approximately  0.22. 

Photoionization  by  solar  EUV  photons  is  the  dominant  production  mechanism  for  O II 
834  A  photons  while  the  photoelectron  impact  ionization  source  contributes  only  about  5%  of  the 
photons.  One  should  note  that  older  values  of  cross  sections  used  in  earlier  studies  [Anderson 
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and  Meier,  1985;  McCoy  et  al.,  1985]  had  given  photoelectrons  a  more  significant  role.  The 
present  calculations  include  the  source  due  to  production  of  834  A  photons  by  photoelectron 
impact  ionization,  since  we  use  the  PEGFAC  code,  which  automatically  incorporates  that  effect. 
Another  minor  contributor  to  the  photon  source  is  the  resonant  scattering  of  solar  834  A 
emissions;  for  our  limb-viewing  geometry,  this  is  a  small  effect  (a  few  percent  of  the  volume 
emission  rate,  j0 ,  for  altitudes  near  zmax)  [Meier,  1991].  As  in  the  tests  of  McCoy  et  al.  [1985] 
and  Anderson  and  Meier  [1985],  our  forward  model  does  not  include  the  effects  of  resonant 
scattering  of  sunlight.  Actually,  the  initial  source  due  to  solar  834  A  emissions  closely  tracks  the 
0+  number  density  distribution  and  tends  to  enhance  the  retrieval  accuracy.  Omitting  this  factor, 
therefore,  renders  our  estimates  of  performance  somewhat  more  conservative.  Of  course,  for 
optimal  accuracy  in  the  retrieval  of  parameter  values  from  real  data,  the  calculation  will  include 
this  source. 

The  forward  model  specifies  a  discrete  grid  of  altitudes  Zj  on  which  to  solve  an  integral 
equation  for  the  volume  excitation  rate  of  line  k  of  the  triplet  [Anderson  and  Meier,  1985; 
Meier,  1991]: 


jk(z)  =  JOk(z)  +  a  Ok  N0+  (Z) 


1  Jjk(z')H 


r  k  ^  k 


*k 


dz' 


ZO 


(10) 


where  the  first  term  on  the  right-hand  side  (RHS)  is  the  initial  photon  production  term,  provided 
by  PEGFAC.  The  second  term  accounts  for  propagation  of  photons  from  other  altitudes  to  a 
given  altitude  z  and  includes  the  effects  of  multiple  resonant  scattering  by  0+.  In  Equation  (10), 
z0  is  the  altitude  below  which  either  the  initial  photon  source  term  is  negligible  or  extinction  is 
large;  a0^  is  the  resonant  scattering  line-center  cross  section  for  line  k;  Tj-(z)  is  the  vertical 
resonant  scattering  line-center  optical  depth  at  altitude  z  ;  t^(z)  is  the  vertical  pure  absorption 
optical  depth,  corresponding  to  absorption  by  N2  ,  O  ,  and  O2;  and  H  is  the  transport  kernel 
[Strickland  and  Donahue,  1970],  The  kernel  H,  also  known  the  Holstein  probability  function,  is 
proportional  to  the  probability  that  a  photon  will  propagate  from  the  region  (z\  z'+dz')  to  the 
region  (z,  z+dz),  undergoing  resonant  scattering  at  that  point  [Holstein,  1947].  Equation  (10) 
assumes  isotropic,  conservative  scattering  and  complete  frequency  redistribution  in  a  plane 
parallel  atmosphere. 

Given  the  volume  excitation  rate,  we  can  then  perform  an  integration  along  the  line  of 
sight  to  predict  the  column  emission  rate  4nl  of  the  radiation  field  in  Rayleighs,  as  measured  at 
the  satellite  [Anderson  and  Meier,  1985;  Meier,  1991]: 

47tl(r,e)  =  10-6fl  jk(r'(s))Tk(r',r)ds  .  (11) 

J  k 

In  Equation  (1 1),  I  is  the  radiance  in  megaphoton  cm'2  sA  ster1;  T^  is  the  probability  that  a 
photon  traveling  along  the  line  of  site  e  from  r'  will  arrive  at  r,  which  in  this  case  is  the  location 
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of  the  observing  instrument.  The  spectral  line  shape,  assumed  here  to  be  Gaussian,  enters 
Equations  (10)  and  (1 1)  implicitly  through  the  integrals  for  H  and  T,  respectively.  The  model 
parameters  which  define  the  0+  number  density  enter  the  forward  model  explicitly  in  Equation 
(10)  and  implicitly  in  both  equations  through  the  resonant  scattering  optical  depth.  Thus  the 
values  of  the  column  emission  rate  on  the  specified  grid  of  tangent  altitudes  (or  OZAs) 
correspond  to  the  model  data  vector  d  in  Equation  (2).  Similarly,  the  integral  on  the  RHS  of 
Equation  (11)  corresponds  to  the  vector  function  G(m),  unless  we  add  a  parameter  m4  to  scale 
the  forward  model  intensity  estimates.  In  this  latter  case,  for  the  i*  line  of  sight  vector,  e;,  G,  = 
m4  x  RHS  of  Equation  (11).  Section  6.3  previews  our  results  when  we  retrieve  the  parameter  m4 
along  with  the  Chapman  parameters  .  The  solution  of  Equation  (10)  is  the  key  factor  in  the 
integrand  of  Equation  (1 1)  and  therefore  constitutes  part  of  the  forward  model.  The  Appendix 
summarizes  the  numerical  approach  to  solving  Equations  (10)  and  (11). 


3.  Model  System  and  Simulated  Observations  for  Tests 

For  the  purpose  of  determining  the  potential  capabilities  of  the  DIT  inversion  technique, 
we  assume  that  the  airglow  is  observed  by  an  instrument  on  a  satellite  at  an  altitude  of  850  km. 
The  systems  of  interest  (e.g.,  RAIDS)  scan  the  limb  from  10°  to  26.5°  below  the  horizon,  so  that 
the  tangent  altitude  varies  between  740  and  90  km,  respectively.  If  the  spacecraft  were 
stationary,  the  distance  of  the  tangent  point  from  the  satellite  would  vary  from  1200  to  3200  km, 
respectively.  To  compensate  for  this  as  much  as  possible  in  the  observing  scheme,  the  systems 
either  scan  downward,  if  facing  aft  of  the  spacecraft,  or  upward,  if  facing  forward,  in  order  to 
achieve  a  tangent  point  trajectory  which  is  as  nearly  perpendicular  to  the  surface  of  the  Earth  as 
possible.  Our  tests  implicitly  assume  that  the  F  region  is  approximately  uniform  over  the 
latitudinal  coverage  of  a  scan.  Section  6  discusses  this  approximation  in  the  context  of  future 
plans  to  account  for  latitudinal  variation  of  the  contribution  function  in  Equation  (11). 

As  our  baseline  case  for  computing  synthetic  observational  data  (synonymous  with 
"observations"  for  the  remainder  of  the  paper),  we  have  represented  the  "true"  0+  number  density 
distribution  as  a  Chapman  layer  with  [Nmax  ,  zmax,  H]  =  [1 .0  x  106  cm-3,  350  km,  60  km]. 
However,  we  have  run  rather  extensive  tests  with  synthetic  data  covering  wide  ranges  of  values 
for  the  Chapman  parameters  and  will  report  on  them  as  needed  in  the  subsequent  sections.  The 
neutral  atmosphere  is  always  the  same  in  our  tests,  given  by  MSIS-86  [Hedin,  1987]  for  the 
parameters  shown  in  Table  1.  In  the  table,  "SZA"  denotes  the  solar  zenith  angle,  Fjo .7  is  the 
value  of  solar  10.7  cm  flux  on  the  previous  day,  and  (Fjo.7)  is  the  81 -day  average  value  of  the 
flux.  For  these  parameters,  the  exospheric  temperature  is  1200  K  ,  which  we  have  used  in 
computing  our  resonant  scattering  cross  sections.  The  latitude  and  longitude  correspond  to  the 
tangent  point  whose  tangent  altitude  is  nearest  zmax  . 

Table  1.  MSIS  Parameters  for  Tests 


Year 

Day 

UT(s) 

AP 

SZA 

f10.7 

<F10.7> 

Latitude 

Longitude 

1980 

81 

61440 

18 

43.49 

188 

163 

32 

254 

To  represent  noise  in  the  observations,  we  have  converted  the  column  emission  rate  in 
Equation  (11)  from  Rayleighs  (R)  to  counts  using  the  sensitivity  of  1  count  R_1  s'1,  which  is 
representative  for  the  various  instruments  to  be  flown  (RAIDS:  0.6  count  R_1  s_1  [Christensen  et 
al.,  1993];  LORAAS  (ARGOS),  SSULI:  1.6  count  R-1  s'1  [McCoy  et  al„  1994b]).  By 
assumption,  the  observation  sample  interval  is  1  s.  We  have  then  generated  Gaussian  random 
numbers  of  characteristic  width  (la)  equal  to  (counts) 1/2,  converted  the  numbers  to  Rayleighs, 
and  added  the  numbers  to  the  column  emission  rates  derived  from  Equation  (11). 

The  altitude  grid  for  the  calculation  of  volume  excitation  rates  consists  of  75  grid  points 
covering  the  altitude  range  from  90  to  990  km.  In  the  direction  of  increasing  altitude,  the 
locations  of  successive  grid  points  increase  exponentially  by  a  constant  factor.  The  OZA  grid  for 
computing  limb  intensity  profiles  consists  of  30  points  at  a  spacing  of  0.28°,  beginning  at  10.28° 
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below  the  horizontal,  and  continuing  with  60  points  at  a  spacing  of  0.135°  at  the  lower  tangent 
altitudes.  This  corresponds  to  the  SSULI  limb  scan  [McCoy  et  al.,  1994b]. 
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4.  Retrieval  of  Chapman  Parameters 

o 

This  section  describes  the  results  of  extensive  tests  that  we  have  performed  on  O  II  834  A 
DIT  inversion  algorithms.  For  ease  of  presentation,  we  introduce  the  term  "ideal"  forward  model 
to  signify  a  forward  model  which  satisfies  the  following:  (1)  the  model  depends  on  no 
approximate  or  unknown  physical  quantities  other  than  the  model  parameters  defining  the 
retrieved  number  density  profile  and  (2)  the  parameterization  of  the  retrieved  number  density 
profile  matches  the  functional  form  describing  the  set  of  ground  truth  profiles.  By  definition,  the 
values  of  the  remaining  physical  parameters  in  the  "ideal"  forward  model  are  identical  to  the  true 
values  characterizing  the  system  being  measured.  In  the  limit  of  negligible  statistical  noise,  the 
intensity  calculated  with  the  ideal  forward  model  represents  an  optimal  fit  to  the  intensity  data 
when  the  ground  truth  number  density  profile  is  entered  into  the  forward  model.  A  model  that  is 
not  ideal  could,  for  example,  use  an  incorrect  instrument  sensitivity  factor  so  that,  given  a  ground 
truth  number  density  profile,  the  model  would  actually  predict  intensity  values  that  were  low  or 
high  by  a  constant  factor. 

To  determine  an  upper  limit  to  the  potential  of  the  DIT  inversion  techniques  and  to 
provide  a  baseline  for  subsequent  refinement  of  the  DIT  algorithms,  the  tests  presented  in 
Subsections  4. 1-4.3  assume  that  the  forward  model  is,  in  fact,  ideal.  Since  the  Chapman  layer 
representation  of  the  ion  number  density  profile  is  being  used,  the  length  of  the  model  parameter 
vector  m  is  M  =  3.  Section  4.4  discusses  the  uncertainties  in  a  priori  inputs  to  the  forward  model 
in  realistic,  nonideal  cases. 

4.1.  Typical  Results  for  the  Three-Parameter  Model 

Figure  1  shows  typical  results  for  our  tests  of  the  three-parameter  model.  The  results  of 
each  test  appear  in  two  diagrams: 

(1)  the  retrieved  (DIT)  0+  number  density  versus  altitude  (solid  line),  along  with  error  bars 
derived  from  Equation  (9)  for  propagation  of  errors,  and 

(2)  synthetic  intensity  data  (centers  of  short  vertical  lines)  with  horizontal  ±VN  error  bars  and 
the  DIT  fit  (solid  line). 

In  both  cases,  the  "true"  values  of  [Nmax  ,  H]  were  the  baseline  values,  [1 .0  x  106  cm'3 , 60  km]. 
For  panels  la  and  b,  zmax  was  250  km,  and  for  panels  lc  and  d,  zmax  was  350  km.  Each 
calculation  used  a  different  set  of  random  numbers  to  simulate  noise.  In  both  cases,  the  initial 
parameter  estimates  used  to  start  the  DIT  iteration  process  were  20%  above  the  "true"  values. 

Comparing  the  information  in  the  plots  and  the  figure  caption  for  the  two  cases,  we  see 
that  the  uncertainties  in  the  retrieved  values  of  Nmax  and  zmax  are  significantly  larger  when  the 
true  zmax  is  250  km  and  that  the  uncertainties  drop  to  minimal  values  when  zmax  increases  to 
350  km.  On  the  other  hand,  the  uncertainty  in  the  scale  height  parameter  is  small  in  both  cases 
(see  figure  caption).  Section  4.3  interprets  these  observations  physically  and  reveals  their 
respective  levels  of  generality  by  studying  the  uncertainties  in  the  retrieved  parameters  as 
functions  of  the  ground  truth  parameter  values.  The  first  figure,  however,  does  support  the 
assertion  that  significant  information  on  the  Chapman  parameters  can  be  retrieved  from  O  II  834 
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A  limb  intensity  measurements,  at  least  when  the  ion  distribution  in  the  F  region  follows 
Equation  (1). 

4.2.  Convergence  Properties 

The  nonlinearity  of  the  intensity  calculation  in  the  forward  model  naturally  requires  that 
the  solution  space  be  searched  for  the  optimal  0+  number  density  profile.  The  simplest  and  often 
the  most  cost-effective  search  procedure  involves  iteration  from  an  initial  estimate  of  the  model 
parameters.  We  have  adopted  this  approach  and  have  thoroughly  tested  the  sensitivity  of  the 
retrieved  parameters  to  the  initial  estimates  of  their  values.  We  began  by  creating  synthetic  data 
for  the  baseline  model  [1  x  106  cm'3,  350  km,  60  km]  with  Gaussian  noise  superimposed.  We 
found  that  the  same  values  of  the  retrieved  parameters  were  obtained  independent  of  the  initial 
estimates.  We  have  duplicated  these  tests  for  a  large  number  of  synthetic  data  sets  based  on  a 
wide  range  of  ground  truth  Chapman  parameters  (ranges  shown  in  Figure  3).  In  all  cases,  we 
have  obtained  results  that  were  independent  of  the  initial  estimates  used  to  start  the  DIT 
calculation.  This  satisfies  the  important  condition  of  robustness  against  choice  of  initial 
parameters.  For  a  given  limb  scan,  the  method  also  converges  to  a  "unique"  maximum  likelihood 
solution.  Here  the  term  "unique"  means  that  the  retrieved  parameters  have  identical  significant 
figures  and  magnitude. 

The  next  test  involves  using  the  same  baseline  parameters  for  the  synthetic  data,  but 
superimposing  different  Gaussian  noise  on  the  synthetic  data  for  each  initial  estimate  of  the 
parameters.  This  time  the  retrieved  parameters  are  different  for  each  initial  guess,  due  to  the 
dependence  on  the  Gaussian  noise,  i.e.,  the  maximum  likelihood  solution  changes.  The  results 
are  summarized  in  Table  2,  which  lists  four  initial  estimates  of  the  parameters  and  the  percentage 
differences  of  the  retrieved  parameter  values  from  the  "true"  values.  Beneath  the  table  appear  the 
"true"  parameter  values  used  to  generate  the  synthetic  data  and  the  maximum  la  uncertainties 
obtained  from  the  covariance  matrix.  Ordinarily,  to  compute  the  approximate  la  values,  one 
would  have  to  perform  these  tests  with  a  large  number  of  synthetic  limb  scans,  each  with 
different  simulated  noise.  The  standard  deviations  of  the  retrieved  parameter  values  would  then 
follow  from  the  averages  of  the  squares  of  the  respective  I  True  -  DIT  I  values.  Instead  we  have 
used  Equation  (7)  for  the  covariance  matrix,  [cov  me]. 
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Table  2.  Convergence  Test  Results  for  Three  Parameter  Model* 


Case 

Initial  Guess  + 

Convergence 

steps 

ITrue  -  DITI,  % 

N 

1  ^max 

zmax 

H 

N 

^max 

zmax 

H 

1 

10.0 

600 

120 

27 

4.4 

0.8 

1.1 

2 

2.0 

600 

120 

8 

2.3 

0.9 

0.0 

3 

1.5 

525 

90 

7 

3.7 

1.8 

0.2 

4 

0.5 

250 

30 

6 

0.9 

0.2 

0.1 

*  The  ground  truth  values  [  Nmax,  zmax,  H]  are  [  10  6  cnr3,  350  km  ,  60  km],  and  maximum  la 

values  are  [  6.0  %,  2.7  %,  1.7  %).  Different  noise  was  used  for  each  run. 

+  The  units  for  Nmax  are  10  6  cm'3,  and  the  units  for  both  zmax  and  H  are  km. 


Figure  2  corresponds  to  initial  estimates  1  and  4  in  Table  2.  The  shape  of  the  initial 
estimates  for  the  intensity  profile  are  quite  different  for  the  two  cases,  as  expected  from  the  work 
of  McCoy  et  al.  [1985].  At  convergence,  the  retrieved  intensity  and  ion  number  density  profiles 
are  in  excellent  agreement  with  the  true  profiles  in  all  cases.  The  1  a  uncertainties  are  consistent 
with  that  picture.  We  have  duplicated  the  above  tests  for  a  large  number  of  synthetic  data  sets 
based  on  a  wide  range  of  ground  truth  Chapman  parameters.  In  all  cases,  we  have  obtained 
results  similar  to  those  in  Table  2  and  Figure  2.  Thus  we  conclude  that  convergence  is  excellent 
and  consistent  with  the  expected  uncertainties  in  the  DIT  result. 

While  our  results  are  quite  encouraging,  we  reiterate  that  we  are  treating  an  ideal  situation 
in  which  the  same  forward  intensity  model  and  the  same  functional  form  for  the  0+  number 
density  distribution  are  used  to  generate  both  synthetic  data  and  DIT  estimates.  The  addition  of 
random  numbers  to  simulate  noise  in  the  data  tends  to  mitigate  this  somewhat,  since  the  smooth 
forward  model  estimates  cannot  match  the  synthetic  data  exactly.  In  a  future  paper,  we  will 
explore  the  potentially  more  difficult  situation  in  which  the  representation  of  the  ion  number 
density  distribution  used  to  generate  synthetic  data  is  different  from  that  used  in  the  DIT  forward 
model. 

4.3.  Uncertainties  in  Retrieved  Model  Parameters 

As  indicated  earlier,  the  standard  deviation  o,  =  ([cov  me]i;)1//2  provides  a  useful  measure 
of  the  uncertainties  in  the  ith  retrieved  model  parameter.  To  generate  approximate  values  of  Oj , 
we  compute  the  forward  model  intensity  gradient  vector,  VG  ,  at  the  location  of  the  "true" 
solution  in  model  parameter  space.  This  gradient  vector  consists  of  the  usual  one-sided 
numerical  derivatives.  We  then  use  VG  to  evaluate  Equation  (7)  for  [cov  me]  ,  of  which  o,2  are 
the  diagonal  elements.  We  have  performed  this  procedure  over  a  three-dimensional  grid  of  "true" 
Chapman  parameter  values.  Note  that  the  results  described  below  apply  only  to  limb  scans  from 
an  altitude  of  850  km  and  covering  OZA  between  100°  and  1 16.5°. 
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Figure  3  shows  two-dimensional  contour  plots  depicting  the  dependence  of  each  Gj  on  the 
"true"  Chapman  parameters  underlying  the  observations.  We  have  selected  only  a  subset  of  two- 
dimensional  planes  within  the  three-dimensional  grid  of  "true"  values  ("ionospheric  states"). 

Over  each  plane,  one  Chapman  parameter  remains  fixed  while  the  others  vary.  We  have  chosen 
to  fix  Nmax  at  1  x  106cnr3  and  H  at  60  km,  both  baseline  values,  while  the  fixed  value  for  zmax  is 
250  km,  chosen  because  of  the  interesting  structure  near  that  value.  Higher  values  of  "true  zmax 
show  parameter  uncertainties  with  progressively  less  structure  and  decreasing  uncertainty  values. 
As  an  example  of  the  use  of  Figure  3,  the  percentage  uncertainty  in  the  retrieved  Nmax  for  the 
baseline  "true"  ionospheric  state  appears  in  Figure  3(a)  or  Figure  3(c)  as  approximately  5%. 

The  behavior  of  the  uncertainty  in  the  retrieved  Chapman  parameters,  as  seen  in 
Figure  3,  follows  from  Equations  (1),  (10),  and  (1 1).  Equation  (1 1)  shows  that  the  behavior  is 
the  result  of  the  interplay  among  the  line-of-sight  path  length,  the  enhancement  of  the  volume 
excitation  rate  as  resonant  scattering  increases  with  the  0+  number  density  in  the  F  region,  and 
the  reduction  in  the  834  A  signal  through  self-absorption  and  pure  absorption,  as  represented  by 
the  transmission  function  T  in  Equation  (1 1).  The  dependence  of  the  volume  excitation  rate  and 
transmission  factor  on  0+  number  density  have  opposing  effects  on  the  intensity,  and  local 
extrema  in  the  uncertainty  could  easily  occur  as  the  true  Nmax  varies.  The  value  of  zmax  directly 
determines  the  altitudes  at  which  resonant  scattering  becomes  appreciable  and  strongly  affects 
the  transmission  factor  linking  the  F  region  with  the  satellite.  The  scale  height  affects  the  shape 
of  topside  limb  intensity  profile  and  the  column  density.  The  latter,  in  turn,  affects  the  optical 
depth  entering  T.  The  transmission  function  itself  can  be  a  significant  source  of  information  on 
the  values  of  the  Chapman  parameters. 

The  remainder  of  this  subsection  uses  such  interpretations  of  Equation  (1 1)  to  reveal 
qualitative  aspects  of  the  dependence  of  the  uncertainty  in  retrieved  Chapman  parameters  on  the 
true  values  underlying  the  data.  Because  the  interplay  among  the  factors  in  Equation  (1 1)  is 
extremely  complex  and  because  space  is  limited,  a  discussion  of  structural  details  of  the  contour 
plots  of  uncertainty  is  beyond  the  scope  of  this  paper.  The  most  important  message  of  Figure  3  is 
that  the  uncertainty  is  small  or  moderate  over  a  significant  region  of  parameter  space.  The 
potential  of  the  834  A  DIT  retrieval  algorithms  is,  therefore,  considerable. 

4.3.1 .  Uncertainty  in  Retrieved  Nmax 

The  first  three  contour  plots  in  Figure  3  show  the  behavior  of  a(Nmax).  In  Figure  3a,  for 
which  the  true  values  of  Nmax  and  zmax  vary  and  the  true  value  of  H  is  60  km,  the  uncertainty  in 
retrieved  Nmax  is  largest  for  low  values  of  zmax  ,  dropping  off  quickly  with  increasing  zmax  .  As 
zmax  increases,  the  distance  along  the  line  of  sight  from  the  observer  to  the  tangent  point  at  zmax 
decreases.  As  a  consequence,  the  optical  path  length  to  the  region  of  the  peak  density  drops  off 
rapidly,  permitting  more  information  (i.e.,  photons)  to  reach  the  sensor.  At  lower  values  of  zmax, 
pure  absorption  becomes  appreciable  and  can  mask  information  coming  from  the  F-layer.  Hence, 
the  lower  the  ionosphere,  the  greater  is  the  uncertainty  in  Nmax. 

The  region  of  uncertainty  values  >  5%  in  Figure  3a  extends  to  higher  true  zmax  as  Nmax 
increases.  In  that  region,  the  contours  indicate  a  plateau  with  a  modest  undulation.  On  the  other 
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hand,  the  uncertainty  is  below  5%  over  much  of  the  diagram.  Upon  inspecting  similar  plots  for 
true  H  between  30  and  90  km,  we  have  found  the  same  qualitative  situation  although  the  peak 
uncertainty  is  larger  at  the  two  extremes  than  in  the  40  -  80  km  range.  Also,  as  H  increases,  the 
region  of  uncertainty  values  >5%  moves  toward  higher  zmax  .  When  Nmax  and  H  are  large, 
resonant  scattering  plays  a  significant  role  not  only  in  enhancing  the  volume  excitation  rate  but 
also  in  reducing  the  magnitude  of  the  transmission  function.  Figure  3a  and  its  relatives  at 
different  values  of  H  indicate  that,  as  Nmax  and  H  increase,  the  inhibition  of  photon  transmission 
to  the  observer  is  the  greater  effect,  reducing  information  received  and  increasing  uncertainty  in 
the  retrieved  parameters. 

Figure  3b  shows  the  uncertainty  in  the  retrieved  value  of  Nmax  as  a  function  of  the  true 
values  of  Nmax  and  H  when  the  true  zmax  =  250  km,  at  which  the  uncertainty  is  greatest.  The 
diagram  shows  that,  for  this  low  peak  height,  the  uncertainty  tends  to  be  greater  at  larger  values 
of  H  for  a  given  Nmax  .  For  a  given  value  of  H,  the  uncertainty  peaks  in  the  vicinity  of  Nmax  =  5 
x  105  .  On  the  other  hand,  similar  plots  for  zmax  >  300  km  show  the  trend  of  the  previous 
paragraph:  the  peak  uncertainty  occurs  in  the  upper  right  comer  of  the  [Nmax,  H]  plane  and  the 
peak  value  decreases  with  increasing  zmax.  Figure  3c  and  its  higher  Nmax  relatives  are  consistent 
with  this  picture. 

4.3.2.  Uncertainty  in  Retrieved  zmax 

Figures  3d,  e,  and  f  show  the  plots  of  o(zmax)  corresponding  parametrically  to  those  of 
a(Nmax)  Figures  3a,  b,  and  c,  respectively.  The  morphologies  of  the  a(zmax)  and  o(Nmax) 
contour  diagrams  are  similar.  Figure  3d  shows  variations  with  true  Nmax  and  zmax  for  true  H  = 

60  km.  For  a  given  value  of  H,  the  error  in  zmax  peaks  in  approximately  the  same  region  of  the 
[Nmax>  zmax]  plane  as  that  for  Nmax  ,  although  the  region  of  high  uncertainty  does  not  extend  as 
far  up  the  Nmax  axis.  In  addition,  as  H  varies  from  30  to  90  km,  the  behavior  is  also  similar  to 
that  of  a(Nmax).  Figures  3e  and  3f  confirm  this  picture,  showing  the  same  qualitative  behavior  as 
those  for  o(Nmax).  The  marked  qualitative  similarities  between  zmax  and  Nmax  occur  because 
both  parameters  describe  the  peak  in  0+  number  density  profile,  and  both  relate  monotonically  to 
the  optical  path  length  between  the  detector  and  the  region  of  the  peak  in  the  ion  number  density 
profile.  Plots  of  two-dimensional  confidence  regions  in  [Nmax,  zmax]  space  confirm  that  Nmax 
and  zmax  are  highly  correlated  [Picone  et  al.,  1994,  1995b;  Kelley  et  ah,  1994],  The  most  notable 
difference  in  a(zmax)  and  o(Nmax)  is  that  the  uncertainty  in  zmax  at  a  given  point  is  smaller  than 
that  of  Nmax.  This  is  due  to  the  close  relationship  between  zmax  and  the  altitude  of  the  peak  in 
the  limb  intensity  profile  over  a  wide  range  of  ionospheric  parameter  values. 

4.3.3.  Uncertainty  in  Retrieved  Scale  Height 

The  behavior  of  a(H),  exemplified  by  Figures  3g,  h,  and  i,  is  qualitatively  different  from 
that  of  the  other  two  parameters,  as  shown  by  a  visual  comparison  of  corresponding  diagrams. 
Unlike  zmax  and  Nmax  ,  which  are  related  primarily  to  the  altitude  and  magnitude  of  the  peak 
intensity,  the  scale  height  determines  the  shape  of  the  normalized  0+  number  density  profile  and 
thus,  the  shape  of  the  topside  limb  intensity  profile.  The  contour  diagrams  in  Figure  3  and  their 
relatives  show  the  uncertainty  in  H  to  be  quite  low  over  the  entire  parameter  space.  This  is 
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primarily  because  the  observations  need  not  sample  far  into  the  medium  to  obtain  H.  Within  the 
context  of  these  small  uncertainties,  the  diagrams  indicate  that  the  uncertainty  in  H  varies 
inversely  with  that  of  Nmax  and  zmax  .  For  large  true  Nmax  and  zmax,  the  uncertainty  in  retrieved 
H  will  decrease  if  the  minimum  OZA  of  the  limb  scan  is  reduced  or  if  the  observation  altitude  is 
increased. 

4.4.  Uncertainties  in  Other  Components  of  the  Forward  Model 

Uncertainty  in  the  forward  model  can  arise  from  several  sources;  examples  are  numerical 
procedures,  instrument  calibration,  estimates  of  the  solar  EUV  irradiance,  and  cross  sections, 
including  those  for  photoionization,  photoelectron  impact  ionization,  resonant  scattering  of 
photons  by  ions,  and  absorption  of  photons  by  neutral  species.  The  Appendix  provides  the 
details  of  our  numerical  procedures.  The  physical  quantities  mentioned  as  examples  appear  as  a 
priori  parameters  which  presently  do  not  vary  during  the  retrieval  process.  A  detailed  study  of 
methods  to  compensate  for  errors  in  these  parameters  and  perhaps  to  retrieve  them  from  the 
intensity  data  is  beyond  the  scope  of  this  paper.  We  can,  however,  offer  a  few  observations  on 
this  topic. 

First,  we  can  retrieve  a  data  scaling  model  parameter  m4  to  address  any  errors  which  have 
the  effect  of  scaling  of  the  intensity  by  a  factor  that  is  approximately  uniform  in  altitude.  This  is 
the  case  if  the  assumed  instrument  sensitivity  is  in  error,  or  if  an  approximately  altitude- 
independent  scaling  error  exists  in  the  g-factor  for  photon  production,  or  if  some  combination  of 
the  two  occurs.  Since  photoionization  of  atomic  oxygen  is  dominant,  a  uniform  error  in  the 
magnitude  of  the  associated  cross  section  should  appear  as  an  approximately  uniform  scaling 
error  in  the  834  A  g-factor.  A  subsequent  paper  [Picone  et  al.,  1995a]  deals  with  such  situations; 
Section  6  previews  our  results. 

The  influence  of  the  various  pure  absorption  cross  sections  requires  a  separate 
investigation.  An  inspection  of  Equation  (10),  however,  allows  us  to  determine  the  ability  to  use 
additional  parameters  to  compensate  for  errors  in  the  resonant  scattering  cross  sections.  Notice 
the  product  of  the  line  center  resonant  scattering  cross  section  aok  and  Nmax  preceding  the 
integral.  The  only  other  factor  containing  ook  is  the  resonant  scattering  optical  depth  Tk ,  which 
is  proportional  to  ook  x  Nmax .  This  means  that  any  model  parameter  scaling  cok  can  be  absorbed 
into  the  model  parameter  Nmax  .  No  method,  therefore,  can  retrieve  both  types  of  parameters 
uniquely  from  O II  834  A  intensity  data  alone.  If  we  have  additional  coincident  data,  such  as  F 
region  electron  density  profiles  measured  by  incoherent  scatter  radar,  and  if  we  have  the  data  in 
sufficient  quantity,  we  should  be  able  to  calibrate  the  DIT-retrieved  value  of  Nmax  .  We  could 
then  derive  a  correction  factor  for  (aok),  where  the  brackets  indicate  an  average  over  the  lines  in 
the  triplet,  or  we  could  derive  correction  factors  for  the  individual  ook  ,  if  the  observations 
resolve  the  834  A  multiplet  components.  One  should  note  that  the  resonant  scattering  cross 
sections  for  O  II  834  A  emissions  are  known  to  within  10%  (Ho  and  Henry,  1983). 
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5.  Concentration  Scalar  Model 


The  previous  section  discusses  uncertainties  in  the  retrieved  Chapman  parameters  as 
functions  of  the  "true"  state  of  the  ionosphere  (i.e.,  the  "true"  parameter  values).  The  desired 
uncertainties  are  often  those  of  the  number  density  at  various  altitudes  within  the  F-region,  as 
retrieved  from  a  given  limb  intensity  data  set.  Because  of  resonant  scattering  and  pure 
absorption,  the  densities  at  lower  altitudes  might  not  contribute  significantly  to  limb  intensity 
data  that  correspond  to  lower  tangent  heights.  The  propagation  of  errors  formula  in  Equation  (9) 
assumes  that  the  low  altitude  extrapolation  provided  by  the  Chapman  layer  is  correct. 

At  the  suggestion  of  R.  Thomas  [private  communication,  1993],  Meier  and  Picone 
[1994]  introduced  the  concentration  scalar  model  (CSM)  to  determine  more  directly  which 
portions  of  the  atmosphere  actually  contribute  to  the  information  received  by  a  UV  sensor 
viewing  the  limb  from  a  satellite.  In  the  CSM,  the  model  parameters  are  the  values  of  the  ion 
number  density  at  the  various  altitudes  on  a  grid.  Unlike  the  Chapman  layer  model  of  Equation 
(1),  the  model  parameters  of  the  CSM  are  not  coupled  through  a  constraining  functional  form. 

As  a  consequence,  the  retrieval  process  must  allow  for  the  possibility  that  the  data  are  not 
sensitive  to  densities  at  some  altitudes;  for  example,  this  could  be  caused  by  extinction  at  the 
lower  altitudes.  We  must,  therefore,  input  a  priori  values  of  the  number  density  and  the 
corresponding  uncertainty  (variance)  at  the  various  grid  points.  For  any  grid  points  from  which  a 
negligible  number  of  photons  reach  the  sensor,  the  DIT  code  will  return  the  respective  a  priori 
values  of  the  number  density.  Similarly,  [cov  me],  the  covariance  of  the  retrieved  model 
parameters,  reverts  to  [cov  m],  the  a  priori  covariance,  at  those  grid  points.  At  other  points,  from 
which  appreciable  information  reaches  the  sensor,  the  variances  of  the  retrieved  density  values 
will  be  smaller  than  those  of  the  a  priori  values.  This  improvement  arises  from  the  added 
information  provided  by  the  data,  as  shown  by  the  first  term  oh  the  RHS  of  Equation  (7).  On  the 
other  hand.  Equation  (6)  shows  that,  if  the  data  are  sufficiently  noisy  (large  variances)  or  if  the  a 
priori  parameter  values  are  known  with  sufficient  accuracy  (small  variances),  then  the  data  will 
contribute  little  to  the  "retrieved"  density  values. 

Our  primary  goals  in  using  the  CSM  in  the  present  paper  are  to  determine  those  altitudes 
from  which  appreciable  amounts  of  information  reach  the  sensor  and  to  assess  the  relative 
importance  of  each  altitude  region  as  compared  to  others.  Thus,  we  require  only  approximate, 
relative  values  of  the  variances,  G;2  ,  of  the  model  parameters  (0+  number  densities),  where  the 
index  i  runs  over  the  altitude  grid.  Our  "standard"  or  baseline  case  for  computing  synthetic 
observational  data  is  again  a  Chapman  layer  with  [Nmax  ,  zmax  ,  H]  -  [1.0  x  106  cm-3,  350  km, 
60  km].  The  number  of  parameters  equals  the  number  of  points  on  the  altitude  grid,  i.e.,  M  =  75. 


To  obtain  the  relative  variances,  we  use  the  same  technique  as  in  Section  4.3  .  First,  we 
compute  the  forward  model  intensity  gradient  vector,  VG  ,  at  the  location  of  the  "true"  solution 
in  model  parameter  space.  Here  the  gradient  vector  consists  of  the  usual  one-sided  numerical 
derivatives.  We  then  use  VG  to  evaluate  Equation  (7)  for  [cov  me],  of  which  Oj2  are  the 
diagonal  elements.  For  the  CSM,  Equation  (7)  requires  a  covariance  matrix,  [cov  m],  for  a  priori 
model  parameter  values.  In  this  case,  [cov  m]  merely  provides  an  upper  bound  on  our  relative 
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variances  and  is  somewhat  arbitrary.  For  the  present  tests,  we  have  assigned  an  upper  bound  of 
20%  to  the  uncertainties  of  the  a  priori  model  parameter  values.  Because  the  concentrations  at 
the  various  altitudes  are  independent  of  each  other,  [cov  m]  is  diagonal,  with  element  ij  is  given 
by  8jj  [0.2  N  +  (z,)]2.  The  symbol  8y  is  the  Kronecker  delta. 


Figure  4  shows  typical  plots  of  the  standard  deviation  of  the  retrieved  number  density 
values,  using  [cov  me]  computed  from  the  CSM.  The  individual  panels  correspond  to  different 
values  of  the  "true"  Chapman  parameters  from  which  the  synthetic  data  were  generated. 
Appearing  above  panels  (b)-(f)  are  the  "true"  parameter  values  which  differed  from  the  baseline 
values  when  the  respective  synthetic  limb  intensity  profiles  were  generated.  The  figure  shows 
that,  with  the  exception  of  very  low  Chapman  layers  (true  zmax  =  250  km),  altitudes  somewhat 
below  the  true  zmax  and  higher  contribute  appreciable  information  for  determining  the  optimum 
result.  The  cases  with  higher  and  lower  Nmax  are  quite  similar  to  the  baseline  case  while  an 
extremely  low  value  of  H  shows  stronger  dependence  on  altitudes  near  zmax. 

Also  shown  in  Figure  4  are  the  uncertainties  computed  from  Equation  (9)  for  the  three- 
parameter  forward  model.  The  Chapman  layer  parameterization  permits  more  accurate  retrieval 
of  the  ion  density  profile,  which  we  expect,  since  our  "observational  data"  are  derived  from  a 
Chapman  layer  representation  of  the  ionosphere.  In  addition,  the  small  number  of  Chapman 
parameters  and  the  dependence  of  the  density  at  each  altitude  on  all  three  parameters  mean  that 
the  Chapman  parameters  can  be  evaluated  without  sampling  the  lower  altitudes. 
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6.  Summary  and  Discussion:  Practical  Applications 

Our  purpose  in  this  paper  has  been  to  demonstrate  the  potential  capability  and  structure 
of  algorithms  based  on  discrete  inverse  theory  for  retrieving  the  distribution  of  0+  in  the  dayside 
F  layer  of  the  ionosphere  from  the  O  II  834  A  airglow.  Our  immediate  application  is  to  satellite- 
based  observations  by  a  limb-scanning  EUV  spectrograph.  The  main  question  regarding  this 
problem  concerns  the  ability  to  acquire  sufficient  information  on  the  0+  distribution  from  the  O 
II  834  A  emissions,  which  are  moderately  optically  thick  due  to  resonant  scattering  within  the  F 
layer.  The  optical  path  length  along  the  line  of  sight  can  be  quite  large  when  viewing  the  limb  of 
the  Earth,  preventing  photons  originating  in  the  lower  thermosphere  from  reaching  the  observing 
system.  To  test  the  potential  capability  of  the  method  we  have  used  the  three-parameter 
Chapman  layer  representation  of  the  0+  distribution  and  a  robust  forward  model  for  computing 
the  limb  intensity,  including  the  effects  of  resonant  scattering  by  0+  and  pure  absorption  by  the 
major  neutral  species. 

In  this  paper,  we  have  treated  only  the  "ideal"  case.  In  the  ideal  case,  the  forward  model 
is  "exact,"  containing  all  of  the  physics  necessary  to  represent  0+  distributions  in  the  F  region 
and  to  compute  accurate  limb  intensity  profiles,  excluding  statistical  noise.  Hence,  our  results  for 
the  ideal  case  represent  an  approximate  upper  bound  on  the  capability  that  we  seek  to  develop. 
For  "nonideal"  cases,  the  DIT  forward  model  differs  materially  from  "reality,"  as  represented  in 
this  case  by  the  model  used  to  generate  "observations."  For  example,  the  values  of  key  cross 
sections  used  in  the  forward  model  have  an  associated  measurement  uncertainty.  Alternatively, 
the  analytic  representation  of  the  0+  number  density  distribution  in  the  retrieval  code  might  not 
be  identical  with  the  functional  form  of  the  ground  truth  distribution.  Section  6.3  below  previews 
recent  results  on  our  approach  to  a  particular  nonideal  case. 

Our  tests  have  shown  that  the  three-parameter  forward  model  always  converges  to  a 
unique  maximum  likelihood  solution  from  a  variety  of  initial  estimates  of  the  Chapman 
parameters,  even  when  the  initial  peak  ion  number  density  is  an  order  of  magnitude  too  large  and 
the  other  two  parameters  are  high  by  a  factor  of  two.  We  have  mapped  the  uncertainty  in  each 
parameter  over  the  three-dimensional  space  covered  by  wide  ranges  of  true  parameter  values. 
Clearly  some  ionospheric  states  are  more  accurately  retrievable  than  others. 

The  concentration  scalar  model  has  provided  a  more  exact  picture  of  the  region  of  the 
ionosphere  being  probed  by  the  limb-scanning  measurements.  The  CSM  shows  that  information 
on  the  ion  number  density  profile  is  available  over  the  entire  top  side  down  to  the  peak,  and 
usually  somewhat  lower.  This  result  gives  us  confidence  in  the  viability  of  the  DIT  technique  for 
deriving  information  on  the  F  region  of  the  ionosphere  from  834  A  limb  scans.  Although  our 
tests  have  been  specific  to  a  single  set  of  parameters  for  satellite-borne  limb  scanning,  we  expect 
the  results  to  apply  to  a  range  of  satellite  altitudes  and  orbits.  In  fact,  the  basic  elements  of  the 
retrieval  algorithms  also  apply  to  measurements  from  sounding  rockets  [Dymond  et  al.,  1995]. 
Coincident  ionosonde  measurements  have  verified  the  0+  number  density  retrieved  by  Dymond 
et  al.,  indicating  that  the  present  paper  provides  a  strong  foundation  for  practical,  robust  methods. 
In  the  following  subsections,  we  mention  some  of  the  issues  which  will  be  addressed  in  future 
publications. 
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6.1.  Assumptions  Regarding  the  Observations  and  the  Ionosphere 

First  we  consider  factors  related  to  the  design  of  the  observing  system.  Equation  (11) 
shows  that  we  have  not  included  the  finite  field  of  view  of  the  sensor  in  our  intensity  integral, 
theoretically  we  would  have  to  average  over  a  range  of  lines  of  sight.  Our  present  tests  therefore 
include  an  implicit  assumption  of  uniformity  across  the  sensor  field  of  view.  For  the  satellite-  ^ 
borne,  limb-scanning  systems  mentioned  earlier,  the  vertical  field  of  view  is  approximately  0.1 
and  is  small  enough  that  our  results  are  not  severely  impacted  by  this  approximation.  Pointing 
accuracy  will  also  affect  the  accuracy  of  the  retrieval  process.  If  the  look  angle  is  off  by  0.1°,  the 
minimum  tangent  height  could  be  off  by  approximately  5  km.  The  retrieved  altitude  of  the  peak 
0+  density  could  be  off  by  a  similar  amount.  The  formalism  of  Section  2  for  handling 
inaccuracies  in  the  forward  model  will  allow  us  to  investigate  this  aspect. 

We  have  treated  the  distributions  of  ions  and  neutrals  as  uniform  over  the  latitude  range 
covered  during  a  scan.  For  the  834  A  triplet,  which  can  be  moderately  optically  thick,  this  is  a 
better  approximation  than  for  an  optically  thin  emission  because  a  significant  amount  of 
information  propagating  from  the  distant  latitudes  can  be  scattered  out  of  the  line  of  sight.  On 
the  other  hand,  this  aspect  should  be  straightforward  to  handle  in  the  DIT  framework  by 
including  a  latitudinal  parameterization. 

Presently  we  evaluate  the  distributions  of  the  major  neutral  constituents  using  the  MSIS- 
86  model  [Hedin,  1987].  In  practice,  for  RAIDS,  ARGOS,  and  SSULI,  we  could  use  the  results 
of  the  DIT  algorithm  for  extracting  such  information  from  FUV  limb  intensity  measurements 
[Meier  and  Picone,  1994].  More  generally  we  might  combine  the  EUV  and  FUV  data  sets  to 
extract  ion  and  neutral  number  density  profiles  simultaneously,  although  the  practical  advantages 
of  such  an  approach  have  not  been  analyzed. 

Anderson  et  al.  [1987]  have  indicated  that  the  simple  Chapman  layer  used  here  is  not  a 
sufficiently  robust  parameterization  of  the  ion  distribution  produced  by  realistic  numerical 
simulations  and  perhaps  the  actual  F  region.  D.  N.  Anderson  and  M.  Fox  [private 
communication,  1993]  have  instead  proposed  more  general  forms  in  which  the  model  0+  profile 
contains  "topside"  and  "bottom-side"  parameters;  the  scale  heights  might  also  be  linearly 
dependent  on  altitude.  The  total  number  of  parameters  would  then  vary  between  M  =  5  and  7. 
We  are  currently  testing  such  representations  against  synthetic  data  based  on  empirical  databases 
and  models. 

6.2.  Numerical  Procedures  and  Solution  Criteria 

A  practical  algorithm  must  consistently  converge  to  acceptably  accurate  representations 
of  the  0+  distribution  in  the  daytime  ionosphere  and  must  do  so  in  an  acceptably  short 
computational  time  interval.  We  have  not  optimized  the  present  algorithms  and  numerical 
procedures  for  high-speed  processing.  The  factors  to  be  considered  are  as  follows:  the  accuracy 
of  the  initial  estimate  of  parameter  values,  definitions  of  convergence  criteria  for  the  iterative 
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process,  time  constraints  placed  on  the  retrieval  activity  by  the  user  and  the  observing  system, 
and  computer  hardware. 

The  best  method  of  increasing  retrieval  speed  is  to  obtain  a  sufficiently  accurate  initial 
estimate  of  the  parameter  vector.  An  example  of  a  strategy  for  achieving  this  is  to  use  the  results 
of  previous  limb  scans  to  initialize  the  retrieval  from  a  given  scan.  This  is  a  sequential 
estimation  technique  proposed  by  Rodgers  [1976]  and  investigated  by  Ortland  et  al.  [1993].  We 
intend  to  study  variations  on  this  theme  and  to  develop  appropriate  enhancements. 

Regarding  the  termination  point  of  the  iterative  process,  we  have  been  computing  a 
maximum  likelihood  solution,  using  a  generalization  of  the  Levenberg-Marquardt  procedure 
described  in  the  Numerical  Recipes  series  [Press  et  al.,  1989].  Currently  the  process  stops  when 
the  change  in  y}  is  smaller  than  a  preset  fractional  value  (1  O'4).  We  have  investigated  other 
termination  criteria  but  have  not  determined  an  optimal  strategy.  In  addition,  the  LM  method  of 
selecting  the  optimal  step  size  for  each  iteration  can  be  computationally  expensive,  due  to  the 
calculation  of  numerical  derivatives  for  the  834  A  inversion. 

The  optimization  of  computer  codes  also  involves  consideration  of  the  requirements  of 
the  user  and  the  user's  particular  computer  host.  For  example,  if  a  scan  of  the  Earth's  limb 
requires  100  s  and  if  each  day  side  scan  is  to  be  processed  for  a  timely  assessment  of  "ionospheric 
weather,"  the  reformatting  and  calibration  of  raw  data,  retrieval  of  0+  distributions,  production  of 
routine  data  products,  and  delivery  to  the  user  must  require  less  than  200  s  per  scan.  (The  200  s 
time  interval  derives  from  the  fact  that  only  dayside  retrievals  are  being  considered.)  Otherwise 
the  data  processing  activity  will  fall  behind  the  rate  at  which  new  data  arrive.  In  addition  to  the 
computational  complexity  of  the  algorithms,  the  speed  of  computer  hardware  and  the 
optimization  of  computer  codes  that  implement  the  algorithms  are  critical  elements  in  achieving 
an  adequate  duty  cycle. 

Most  importantly,  better  criteria  might  exist  for  defining  the  DIT  solution  than  that  of  the 
maximum  likelihood.  For  the  retrieval  of  a  three-parameter  Chapman  layer  from  a  given  limb 
intensity  profile,  the  DIT  procedure  converges  arbitrarily  close  to  a  unique  maximum  likelihood 
solution.  For  some  ionospheric  states,  however,  the  maximum  likelihood  solution  itself  can 
change  by  a  significant  amount  as  the  particular  realization  of  statistical  noise  changes.  In 
preliminary  tests,  we  have  had  some  success  in  finding  a  different  criterion  for  defining  the  DIT 
solution  so  that  the  solution  varies  less  strongly  with  the  particular  sample  of  noise.  This  issue 
merits  a  careful  and  detailed  study. 

6.3.  Advanced  Topics 

We  conclude  with  a  preview  of  results  on  two  topics  that  extend  beyond  the  present 
baseline  study.  While  we  will  be  reporting  in  detail  on  these  topics  in  succeeding  publications 
[Picone  et  al.,  1995a, b],  interest  has  been  sufficiently  high  to  warrant  mention  here. 
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6.3.1.  Different  Chapman  Layers  With  "Nearly  Identical"  Intensity  Profiles 

Link  et  al.  (1994)  have  computed  model  834  A  limb  intensity  profiles  for  a  satellite  in  an 
850  km  orbit  while  viewing  (model)  ionospheres  with  values  of  Nmax  differing  by  a  factor  of  2 
and  zmax  differing  by  50  km.  Excluding  the  effects  of  noise,  such  limb  intensity  profiles  can 
differ  by  approximately  3%  or  less  at  look  angles  corresponding  to  tangent  altitudes  below  650 
km.  Given  that  this  difference  is  much  less  than  the  scatter  seen  in  data  from  the  STP78-1 
mission  [McCoy  et  al.,  1985],  Link  et  al.  conclude  that,  for  STP78-1  data,  the  uncertainties  in 
Nmax  and  zmax  for  the  retrieved  0+  number  density  profiles  can  be  worse  than  the  respective 
factors  of  2  and  50  km  by  which  the  model  ionospheres  differ. 

Presently  we  are  studying  the  statistical  uncertainty  in  the  retrieved  parameters  in  the 
context  of  future  EUV  remote  sensing  missions,  for  which  the  instruments  are  more  sensitive 
than  that  of  STP78-1  by  a  factor  of  three  (RAIDS)  or  more  (SSULI,  ARGOS).  In  contrast  to  the 
procedure  of  Link  et  al.,  which  was  based  on  inspection,  the  DIT  approach  provides  an  analytic 
method,  Equation  (7),  for  estimating  the  uncertainties  in  the  retrieved  profiles.  In  addition,  the 
minimization  of  a  characteristic  function  (like  y^)  provides  the  basis  for  selecting  an  optimal  set 
of  model  parameters,  thereby  distinguishing  among  various  "nearly  identical"  model  intensity 
profiles  and  their  underlying  0+  number  density  profiles.  This  discrimination  is  statistical 
because,  for  a  given  "true"  ionospheric  state,  the  retrieved  parameter  values  vary  with  the 
particular  realization  of  statistical  noise  present  in  a  given  limb  intensity  profile.  The  statistical 
spreads  in  retrieved  parameter  values,  therefore,  depend  on  the  sensitivity  of  the  observing 
instrument. 

To  verify  this  capability,  we  have  studied  two  different  cases  which  are  similar  to  the 
single  case  shown  by  Link  et  al.  Each  test  case  consists  of  two  "true"  ionospheric  states  with 
significantly  different  values  of  Nmax  and  zmax  but  with  "nearly  identical"  intensity  profiles, 
excluding  statistical  noise.  To  each  intensity  profile  in  a  pair,  we  added  twenty  different  sets  of 
random  noise,  generating  forty  synthetic  limb  intensity  observations  per  pair  of  intensity  profiles; 
we  did  this  for  each  of  three  instrument  sensitivities:  0.2,  0.6,  and  1.0  counts/R/s.  We  then 
retrieved  the  maximum  likelihood  solution  for  each  limb  observation.  Our  results  are  as  follows: 

(1)  Compared  to  the  case  of  low  sensitivity  instruments  like  that  on  STP78-1,  the  uncertainties 
in  the  retrieved  parameters  are  much  lower  for  the  more  sensitive  instruments  aboard  future 
missions.  The  scatter  in  the  noise  is  much  lower,  and  hence,  the  ability  to  resolve  "nearly 
identical"  limb  intensity  profiles  is  significantly  greater. 

(2)  The  cases  in  which  significantly  different  Chapman  layers  yield  nearly  identical  limb 
intensity  profiles  are  those  for  which  the  uncertainty  in  the  retrieved  Nmax  and  zmax  is  high. 

Over  a  great  range  of  values  for  the  "true"  Chapman  parameters,  this  statistical  uncertainty  is  less 
than  5%  (Figure  3),  so  that  the  issue  of  nearly  identical  limb  intensity  profiles  is  much  less 
important. 

(3)  For  future  missions,  the  68%  confidence  regions  of  the  retrieved  solutions  are  well  separated 
and  are  readily  distinguishable. 
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Thus,  we  find  from  this  preliminary  analysis  that,  by  employing  powerful  inversion  methods, 

o 

such  as  DIT,  remote  sensing  of  the  dayside  ionosphere  using  O  II  834  A  measurements  can 
indeed  yield  reliable  0+  concentrations  with  relatively  low  uncertainties. 

6.3.2.  Retrieval  of  a  Data  Scaling  Factor 

For  a  3-parameter  representation  of  the  0+  number  density  profile,  we  have  also 
investigated  extensively  the  properties  of  an  algorithm  which  retrieves  an  additional  parameter, 
1TL4,  that  scales  the  intensity  estimate  by  the  forward  model  prior  to  comparison  with 
observational  data  [Picone  et  al.,  1995a].  This  approach  can  mitigate  the  influence  of  inaccurate 
instrument  calibration,  which  may  be  attributable  to  errors  prior  to  launch  or  to  changes  that 
occur  over  the  course  of  the  mission.  Such  inaccuracy  can  cause  incorrect  scaling  of  the 
intensity,  as  estimated  by  the  forward  model.  Another  factor  which  might  also  be  mitigated  in 
this  way  is  an  error  in  the  magnitudes  of  the  g-factors  (excitation  rate  per  second  per  ion)  for 
production  of  the  O II  834  A  emissions.  Errors  in  g-factors  can  arise  from  uncertainties  in  the 
excitation  cross  sections  and  in  the  solar  flux,  the  neutral  atmosphere,  and  the  photoelectron  flux. 
The  retrieval  of  a  single  data  scaling  parameter  would  be  most  effective  when  the  g-factor  is  off 
by  a  multiplicative  scalar  that  is  approximately  constant  in  altitude  over  the  range  of  interest. 

Our  preliminary  tests  consider  Chapman  layers  with  true  [Nmax,  H]  =  [2.0  x  10^  cm'3, 

60  km]  and  true  zmax  e  {250,  350,  450,  550}  km.  To  obtain  synthetic  observational  data  for  our 
tests,  we  have  generated  intensity  values  with  the  forward  model  and  then,  prior  to  adding 
simulated  noise,  we  have  scaled  them  by  a  factor  which  differed  from  1 .0  (either  1.3  or  0.7).  For 
each  value  of  the  true  zmax,  we  have  used  twenty  different  synthetic  limb  intensity  profiles  (i.e., 
twenty  sets  of  random  noise).  Without  the  retrieval  of  the  fourth  model  parameter,  the  DIT  code 
could  not  obtain  useful  estimates  of  the  true  Chapman  parameters.  However,  with  the  additional 
retrieved  parameter,  the  performance  is  consistent  with  the  previous  sections,  which  have 
assumed  the  ideal  case  (true  1114  =  1 .0).  For  zmax  >  450  km,  the  average  values  of  the  retrieved 
Chapman  parameters  and  the  scaling  factor  are  within  a  few  percent  of  the  true  values,  and  the 
standard  deviations  are  small.  As  the  value  of  zmax  decreases  below  450  km,  the  uncertainty 
increases  as  in  the  ideal  cases  treated  earlier  (Figure  3).  Since  the  true  value  of  014  is  independent 
of  the  state  of  the  local  ionosphere,  one  has  the  luxury  of  selecting  only  limb  intensity  data  which 
satisfy  stringent  requirements.  Note  that  the  systems  considered  here  provide  approximately  30 
limb  intensity  profiles  per  orbit  —  more  than  enough  for  this  purpose.  Once  a  value  for  m4  has 
been  obtained,  one  can  perform  routine  3-parameter  inversions  with  the  forward  model  intensity 
scaled  by  that  value  of  m4.  Our  expectation  is  that  this  technique  will  work  well  for  generalized 
representations  of  the  ionosphere,  for  which  M  >  3.  Since  the  calibration  of  the  instrument  might 
change  over  time,  periodic  updates  of  the  value  would  be  desirable. 
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Appendix:  Numerical  Details  of  the  Forward  Model 

We  briefly  describe  the  forward  model  used  for  these  tests.  The  model  includes  the 
effects  of  multiple  resonant  scattering  by  0+  and  absorption  of  photons  by  O,  N2,  and  O2  .  The 
integral  equation  (10)  is  identical  to  that  discussed  in  the  original  paper  on  the  topic  by  Anderson 
and  Meier  [1985].  Earlier,  Strickland  and  Donahue  [1970]  and  Strickland  and  Rees  [1974]  had 
analyzed  the  general  problem  of  multiple  resonant  scattering  under  the  assumption  of  complete 
frequency  distribution.  Anderson  and  Meier  [1985]  later  modified  the  approach  for  834  A 
airglow,  adding  an  approximate  "external  source"  (boundary)  term  to  Equation  (10)  to  account 
for  the  fact  that  the  region  in  which  the  emissions  are  produced  can  lie  well  below  the  F  layer. 
When  comparing  recent  calculations  to  those  in  the  original  work,  one  must  remember  that  the 
accepted  values  of  cross  sections  for  photon  production,  resonant  scattering,  and  absorption  have 
changed  as  better  data  have  become  available  [Link  et  al.,  1994]. 

The  need  for  the  external  source  term  in  the  method  of  Anderson  and  Meier  arose  because 
the  problem  is  solved  on  the  resonant  scattering  optical  depth  grid,  which  ceases  significant 
variation  with  decreasing  altitude  once  the  "bottom"  of  the  0+  layer  is  reached.  This  requires  the 
establishment  of  an  arbitrary  definition  for  the  lower  boundary  of  the  0+  layer.  Since  most  of  the 
photons  are  usually  produced  below  the  0+  layer,  the  external  source  must  faithfully  represent 
propagation  of  the  photons  into  the  F  region  from  below.  The  code  used  for  the  present  work 
solves  Equation  (10)  on  an  altitude  grid  that  extends  below  the  region  where  834  A  photons  are 
produced  so  that  no  boundary  term  is  required.  For  this  reason,  no  interpolation  of  physical 
quantities  onto  an  optical  depth  grid  is  necessary,  and  the  only  other  approximations  are  purely 
numerical,  such  as,  the  treatment  of  physical  quantities  as  discrete. 

The  user  can  select  the  altitude  grid,  denoted  by  Zj  for  i  ranging  from  1  to  N.  We  are 
presently  using  a  grid  of  75  points,  the  values  of  which  increase  exponentially  from  cell  to  cell; 
tests  have  shown  only  weak  dependence  on  the  particular  gridding  scheme.  We  evaluate  the 
volume  excitation  rates  at  these  altitude  values.  In  the  current  implementation,  we  treat  the  set 
{Zj ;  i  =  1,2, ... ,  N}  as  the  cell  boundaries,  with  cell  centers  at  {zcj  =  xh{zx  +  zl+\),  i  =  1,  2, ... , 
N-l  }.  In  the  following  equations,  the  subscript  i  with  no  superscript  will  label  a  quantity  which 
is  evaluated  at  zx ,  and  the  subscript  i  with  superscript  c  will  label  a  quantity  which  is  evaluated  at 
zcj .  We  then  solve  the  following  numerical  approximation  for  Equation  (10)  at  cell  centers: 


ji  =  joi  (z)  +  °oWo+)cx 
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(A2) 


,  M  J_  7  e-2x2E2(ae-x2+b)  _dx 

£(a>b)  =  2Vi  Ia0No.(zf)e-x2+  £oaNa(zf) 

a 

The  second  term  on  the  RHS  of  Equations  (10)  and  (Al)  includes  multiple  scattering 
effects.  In  Equation  (A2),  the  index  a  ranges  over  the  major  neutral  species,  which  absorb  the 
834  A  triplet  emissions.  The  functional  E2  is  the  second  order  exponential  integral  [Abramowitz 
and  Stegun,  1972].  To  compute  e(a,b),  we  have  used  the  Gauss-Hermite  quadrature  formula  with 
terms  up  to  n  =  20  [Press,  et  al.,  1989;  Abramowitz  and  Stegun,  1972;  Krylov,  1959]. 

For  the  purpose  of  testing  the  DIT  approach,  the  numerical  approximations  have  been 
straightforward.  The  following  assumptions  are  implicit  in  the  above  equations: 

(1 )  the  cell  sizes  are  small  enough  so  that  the  volume  excitation  rate  is  approximately 
constant  over  each  cell; 

(2)  the  variation  with  altitude  of  both  the  resonant  scattering  and  pure  absorption  optical 
depths  is  approximately  linear  about  the  center  of  each  cell; 

(3)  the  lower  bound  of  all  frequency  integrations  for  the  Gaussian  line  shape  can  be 
approximated  by  -  <»  ,  incurring  negligible  errors. 


To  solve  Equation  (Al),  we  subtract  the  second  term  on  the  RHS  from  both  sides,  giving  us  a 
matrix  equation,  [A]  j  =  j0 .  One  can  then  use  a  standard  matrix  inversion  technique  to  solve  for 
the  total  volume  excitation  rate  j,  which  rapidly  approaches  j0  below  the  Chapman  layer.  At  each 
look  angle,  we  compute  the  intensity  in  Equation  (1 1)  by  using  a  Simpson's  rule  integration  with 
60  points  along  the  line  of  sight,  assuming  a  spherical  Earth. 

Comparisons  with  the  Monte  Carlo  model  of  Meier  and  Lee  [1982],  modified  for  the  0+ 
834  A  problem,  have  shown  good  agreement  for  total  volume  excitation  rates  over  the  entire 
range  of  Chapman  parameter  space  covered  by  the  tests.  The  PEGFAC  code  of  Stockland  and 
Meier  [1982],  with  updated  cross  sections  [Meier,  1991;  Link,  private  communication,  1 994] , 
provides  the  initial  source  term,  j0 ,  in  Equation  (Al).  We  are  presently  developing  a  synthesis  of 
the  initial  source  term  (i.e.,  the  g-factor)  based  on  recent  theoretical  work  by  Meier  [Majeed  et 
al.,  1993]. 
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Figures 


1 .  Results  of  DIT  retrieval  for  a  Chapman  layer  with  true  [Nmax  ,  H]  =  [1 .0  x  106  cm’3,  60 
km]  and  zmax  =  (a),  (b)  250  km  and  (c),  (d)  350  km.  Panels  (a)  and  (c)  show  respective 
profiles  of  the  DIT  retrieved  0+  number  density  versus  altitude  (solid  line)  and  error  bars 
derived  from  Equation  (9).  The  retrieved  parameters  and  lo  uncertainties  are  [0.94  x  106 
cm*3,  255  km,  61  km]  and  [33%,  2.9%,  1.9%]  for  zmax  =  250  km  and  [0.98  x  10*  cnr3, 
356  km,  60  km]  and  [4.1%,  2.4%,  1.8%)  for  zmax  =  350  km.  Panels  (b)  and  (d)  show 
respective  synthetic  intensity  data  (short  vertical  lines  with  horizontal  ±  Vn  error  bars) 
and  the  DIT  intensity  fit  (solid  line). 

2.  Retrieved  0+  number  density  profiles  and  corresponding  intensity  profiles  for  cases  1  and 
4  of  Table  2  with  true  parameters  given  by  the  baseline  case  [Nmax  ,  zmax  ,  H]  =  [1 .0  x 
106  cm-3,  350  km,  60  km].  Diagrams  (a)  and  (b)  depict  typical  results  when  starting  the 
iterative  DIT  process  with  Chapman  parameters  that  are  "high"  [1.0  x  107  cnr3, 600  km, 
120  km].  Diagrams  (c)  and  (d)  correspond  to  "low"  initial  estimates  [5.0  x  105  cnr3, 250 
km,  30  km].  The  dashed  lines  show  the  initial  DIT  ion  concentration  and  intensity 
profiles;  the  solid  lines  show  the  final  DIT  results;  the  dotted  lines  show  the  true  ion 
density  profiles;  and  the  plus  (+)  signs  show  the  synthetic  data,  generated  with  different 
random  noise  for  the  two  cases. 

3.  Contour  plots  of  the  la  uncertainties  in  the  retrieved  [Nmax  ,  zmax  ,  H],  with  each  la 
value  being  expressed  as  a  percentage  of  the  respective  true  parameter  value.  The  first 
three  plots  show  a(Nmax)  as  a  function  of  the  true  values  of  (a)  Nmax  and  zmax  when  the 
true  scale  height  H  =  60  km,  (b)  Nmax  and  H  when  the  true  zmax  =  250  km  and  (c)  zmax 
and  H  when  the  true  Nmax  -  1.0  x  106  cnr3.  Plots  (d),  (e),  and  (f)  show  a(zmax)  and 
correspond  to  (a),  (b),  and  (c),  respectively.  Similarly  plots  (g),  (h),  and  (i)  show  a(H). 

4.  Typical  plots  of  the  standard  deviation  of  the  retrieved  number  density  values  for  the 
concentration  scalar  model  as  compared  to  the  three-parameter  Chapman  layer  model, 
when  true  values  are  given  by  (a)  baseline  case,  [Nmax  ,  zmax  ,  H]  =  [1.0  x  106  cnr3,  350 
km,  60  km],  and  (b)  -  (f)  baseline  values  with  the  exception  of  the  parameter  indicated. 
The  +  signs  indicate  CSM  results  for  the  respective  cases.  Dotted  lines  indicate  the 
uncertainties  resulting  from  a  three-parameter  Chapman  layer  representation. 
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Case  1  :  High  Initial  [  N 


Number  Density  (10°  cm  J)  834 


(b)  H  =  30  km 
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